#!/usr/bin/perl -w
use strict;
use warnings;
use FindBin;
use lib ($FindBin::Bin, "$FindBin::Bin/../perl_modules");
use die_error;
use Getopt::Long;

my $usage = <<'USAGE';

USAGE:

parse_blastall_output.pl < blast_all.results
				[--minimum_alignment INTEGER]
				[--reciprocal = 0]
				[--debug]
				[--help]

    Parses blastall results in XML format (alignment view option -m 7)
	into the standard format produced by seq_pairs_align
	--reciprocal parse results even if id1 == id2. Otherwise only id1 < id2 is accepted
	--minimum_alignment is the minimum alignment result which would be accepted.
	Defaults to 50
USAGE

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#   Get options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
# mandatory


# optional parameters
my $help = undef;
my $reciprocal = 1;
my $minimum_alignment = 50;
my $debug = undef;
GetOptions(
			'help'			   		=> \$help,
			'minimum_alignment=s'	=> \$minimum_alignment,
			'reciprocal=i'			=> \$reciprocal,
			'debug=s'				=> \$debug,
			);

die $usage if ($help);


if ($debug)
{
	open_or_die(*DEBUG,  ">$debug");
}

my %data;
my $prot_name = '';
sub do_print($$\%);
my $cnt_hits = 0;
my $cnt_prot = 0;
my ($i_len, $i_name);
while (<>)
{
	chomp;
	if (/^<\/BlastOutput>$/)
	{
		print "#Blastall EOF\n"
	}
	if (/<([^>]+)>([^<]+)</)
	{
		if ($1 eq 'Hsp_midline')
		{
			$cnt_hits++;
			$data{$1} = $2;
			if ($data{'Iteration_query-def'} ne $prot_name)
			{
				$cnt_prot++;
				$prot_name = $data{'Iteration_query-def'};
			}
			$i_len	= $data{'Iteration_query-len'};
			$i_name = $data{'Iteration_query-def'};

			# only save top hit
			if (exists $data{'Hit_def'})
			{
				do_print($i_len,$i_name, %data);
			}
			%data = ();
			$data{'Iteration_query-len'} = $i_len;
			$data{'Iteration_query-def'} = $i_name;
		}
		else
		{
			$data{$1} = $2;
		}
	}
}
#print STDERR "hits = $cnt_hits\n";
#print STDERR "prot = $cnt_prot\n";

sub gap_positions_array_to_str(\@)
{
	return join ("|", map{$_->[0].':'.$_->[1]} @{$_[0]});
}
sub do_print($$\%)
{
	my ($i_len, $i_name, $data) = @_;
	my ($j_len, $j_name) = @$data{'Hit_len', 'Hit_def'};
	($i_name) = split / /, $i_name;
	($j_name) = split / /, $j_name;
	return if (!$reciprocal && $i_name ge $j_name);
	if ($debug)
	{
		print DEBUG "$i_name\t$j_name\n";
		print DEBUG "\t".substr($data->{'Hsp_qseq'}, 0, 100), "\n";
		print DEBUG "\t".substr($data->{'Hsp_hseq'}, 0, 100), "\n";
		print DEBUG "\t".substr($data->{'Hsp_midline'}, 0, 100), "\n";
	}
	my ($i_aa_line, $j_aa_line) = @$data{'Hsp_qseq', 'Hsp_hseq'};
	my (@i_pos,  @j_pos);
	my ($i_beg,  $j_beg, $i_end,  $j_end)
				= @$data{'Hsp_query-from', 'Hsp_hit-from','Hsp_query-to', 'Hsp_hit-to'};
	$i_beg--;
	$j_beg--;


	# note the offset to the first aligned position
	my ($align_beg) = 0;
	my ($align_len) = $data->{'Hsp_align-len'};
	return if ($align_len <= $minimum_alignment);
	if ($i_beg > $j_beg)
	{
		push(@j_pos, [0, $i_beg - $j_beg]);
		$align_beg = $i_beg;
	}
	elsif ($i_beg < $j_beg)
	{
		push(@i_pos, [0, $j_beg - $i_beg]);
		$align_beg = $j_beg;
	}

	# save the number and position of the aligning characters for each sequence
	while ($i_aa_line =~ /(\-+)/g)
	{
		# save position and length of aligning markers
		my $len = length($1);
		my $pos = pos($i_aa_line) - $len + $i_beg;
		push(@i_pos, [$pos, $len]);
		# adjust starting position to account for the markers already counted
		$i_beg -= $len;
	}
	while ($j_aa_line =~ /(\-+)/g)
	{
		# save position and length of aligning markers
		my $len = length($1);
		my $pos = pos($j_aa_line) - $len + $j_beg;
		push(@j_pos, [$pos, $len]);
		$j_beg -= $len;
	}
	my $i_len_final = $i_len;
	$i_len_final += $_->[1] for (@i_pos);
	my $j_len_final = $j_len;
	$j_len_final += $_->[1] for (@j_pos);


	# add end alignment marks
	my $final_len = $i_len_final;
	if ($i_len_final > $j_len_final)
	{
		push(@j_pos, [$j_len, $i_len_final - $j_len_final]);
		$final_len = $i_len_final;
	}
	elsif ($j_len_final > $i_len_final)
	{
		push(@i_pos, [$i_len, $j_len_final - $i_len_final]);
		$final_len = $j_len_final;
	}

	my $i_gap_str = gap_positions_array_to_str(@i_pos);
	my $j_gap_str = gap_positions_array_to_str(@j_pos);
	my $residues = $align_len;
	my @i_aa = split (//,  $i_aa_line);
	my @j_aa = split (//,  $j_aa_line);
	my $cnt_identical = 0;
	for (0..$#i_aa)
	{
		++$cnt_identical if ($i_aa[$_] eq $j_aa[$_]);
	}
	$cnt_identical /= @i_aa;
	$residues -= $i_aa_line =~ tr/\-/-/;
	$residues -= $j_aa_line =~ tr/\-/-/;
	my $max_len = ($i_len > $j_len) ? $i_len : $j_len;
	
#seq1    seq2    start   len     positions       coverage        e_value len1    len2    residues        percent_id      bitscore  gaps1    gaps2

	print join("\t",  $i_name, $j_name, $align_beg, $align_len, 
					$final_len, # positions
					$residues / $max_len, 
					$data->{'Hsp_evalue'}, 
					$i_len, 
					$j_len,
					$residues, 
					$cnt_identical,
					$data->{'Hsp_bit-score'}, 
					$i_gap_str, $j_gap_str), "\n";
}





=pod
Iteration_query-def			ENSAAAP0000002									pep:known									chromosome:NCBI35:1:10370:19200:-1									gene:ENSG00000146556									transcript:ENST00000356572
Hit_def							ENSP00000348979									pep:known									chromosome:NCBI35:1:10370:19200:-1									gene:ENSG00000146556									transcript:ENST00000356572
Hit_len							199 
Hit_num							1   
Hsp_align-len					194 
Hsp_bit-score					340.117
Hsp_evalue						5.05788e-94
Hsp_gaps						26  
Hsp_hit-from					5   
Hsp_hit-to						180 
Hsp_hseq						RRHHDSCEDAALPGRSDLCRALHPARPAARGGRPADGGCPAVPAEGLWRHLQQVGAATHPHLMRGPSLSW------------------AIPATHPLGAAPLGNPGPGEPIWLPCVGGAEVRAQLGSFYSWESESSLPSTPLSLLVSCGGSALTTHAFPPIGGPIPSRFPAANKYLPSSARVQGSVGISIYPATD
Hsp_qseq						RRHHDSCEDAALPGRSDLCRA--------RGGRPADGGCPAVPAEGLWRHLQQVGAATHPHLMRGPSLSWWWWWWWWWWWWWWWWWWWAIPATHPLGAAPLGNPGPGEPIWLPCVGGAEVRAQLGSFYSWESESSLPSTPLSLLVSCGGSALTTHAFPPIGGPIPSRFPAANKYLPSSARVQGSVGISIYPATD
Hsp_query-from					1   
Hsp_query-to					186 
Hsp_score						871 
=cut
